Ground state optimization and hysteretic demagnetization: 
the random-field Ising model 
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We compare the ground state of the random-field Ising model with Gaussian distributed random 
fields, with its non-equilibrium hysteretic counterpart, the demagnetized state. This is a low energy 
state obtained by a sequence of slow magnetic field oscillations with decreasing amplitude. The main 
concern is how optimized the demagnetized state is with respect to the best-possible ground state. 
Exact results for the energy in d = 1 show that in a paramagnet, with finite spin-spin correlations, 
there is a significant difference in the energies if the disorder is not so strong that the states are 
trivially almost alike. We use numerical simulations to better characterize the difference between 
the ground state and the demagnetized state. For d > 3 the random-field Ising model displays a 
disorder induced phase transition between a paramagnetic and a ferromagnetic state. The locations 

(DS) (GS) 

of the critical points Rc \ R c differ for the demagnetized state and ground state. Consequently, 
it is in this regime that the optimization of the demagnetized stat is the worst whereas both deep 
in the paramagnetic regime and in the ferromagnetic one the states resemble each other to a great 
extent. We argue based on the numerics that in d = 3 the scaling at the transition is the same in 
the demagnetized and ground states. This claim is corroborated by the exact solution of the model 
on the Bethe lattice, where the _R c 's are also different. 

PACS numbers: 



I. INTRODUCTION 

The relation between equilibrium and non-equilibrium 
states is a central problem in the physics of disordered 
systems. Disorder induces a multitude of metastable 
states in which the system can easily be trapped. The 
dynamics is usually very slow, or glassy, and on obser- 
vational timescales the system is basically always out of 
equilibrium. On the other hand, from the theoretical 
point of view it is easier to consider equilibrium proper- 
ties, since in this case it is possible to use all the ma- 
chinery of statistical physics to tackle the problem. The 
question is whether the equilibrium properties of disor- 
dered systems provide a faithful representation of the 
non-equilibrium states in which the system is likely to 
be found in practice. This dichotomy is at the core of 
many unsolved issues in the field of disordered system. 
Typical quantities that one could compare are the energy, 
the geometric characterization of the state (as domains 
in magnets), and the energy cost of excitations. 

A simplification of the problem is obtained considering 
only athermal processes, in which the temperature of the 
system plays no role and can be ignored. The equilib- 
rium state is in this case just the ground state (GS), the 
state of minimal energy A zero temperature, non- 
equilibrium dynamics is purely relaxational: the system 
falls simply in the closest metastable state. A convenient 
way to allow the system to explore the various metastable 
states is by applying an external magnetic field. Differ- 
ent field histories typically result in hysteresis and lead 



to different metastable configurations Q. 

The demagnetization process consists in applying a 
slowly varying AC field with decreasing amplitude, and 
provides a simple way to access low energy states 0. 
It has been studied for more than a century, but until 
recently the question how close the demagnetized state 
(DS) is to the true GS was not addressed. This is the 
concern of our work, the problem of how such an opti- 
mization process works in the case of a random magnet. 
Recently, Pazmandy et al. have proposed the demag- 
netization process as the basis for a new optimization 
algorithm for disordered systems The idea behind 
such "hysteretic optimization" , is that demagnetization 
leads to a low energy state, sufficiently close to the GS, 
which can then be reached by applying other methods 
using the DS as an input. The method was tested for 
different models like spin glasses and NP-hard problems. 

Here we will concentrate on the random-field Ising 
model (RFIM), that, while retaining some complex fea- 
tures characteristic of disordered systems, still allows for 
a theoretical analysis Q|. In the RFIM, due to the ab- 
sence of frustration, the equilibrium state is relatively 
simple, however, the non-equilibrium dynamics is far 
from trivial. Due to the coupling of the local disorder 
to the order parameter, even the GS presents a variety of 
phenomena, which can be studied numerically 0, 0,013 • 
In fact the GS of the RFIM can be found in a poly- 
nomial CPU-time, with exact combinatorial algorithms 
and solved exactly in d = 1 and on the Bethe lat- 
tice 0, llfjj . The equilibrium critical exponents for ran- 



2 



dom field magnets have been measured experimentally in 
Fe . 93 Zn . 07 F 2 0G2. 

The hysteretic properties of non equilibrium RF1M 
have been widely studied in the recent literature. The 
hysteresis loops display a disorder induced phase transi- 
tion: for low disorder the loop has a macroscopic jump 
at the coercive field, while at high disorder the loop is 
smooth, at least on the macroscopic scale Q, Il5| . 
At smaller scale the magnetization curve is highly discon- 
tinuous, showing Barkhausen-type bursts, in correspon- 
dence to jumps between different metastable states |l6j . 
A disorder induced non-equilibrium phase transition in 
the hysteresis loop has been studied experimentally in 
Co-CoO films 17| and Cu-Al-Mn alloys 0. 

Extensive numerical simulations have been used to 
characterize disorder induced transitions in the non- 
equilibrium RFIM and critical exponents have been es- 
timated in several dimensions ^j, IT^, l2(ij ]. The model 
has been studied by the renormalization group and the 
exponents have been computed in a e = 6 — d expan- 
sion In addition the hysteresis loop has been com- 
puted exactly in d — 1 and on the Bethe lattice, where 
the disorder induced transition is present for sufficiently 
large coordination number. While in d = 1 there is def- 
initely no transition, the situation in d = 2 is less clear. 
Recently the problem of minor loops has been tackled 
analytically and numerically. In particular, the demag- 
netization curve has been computed exactly in d = 1 [2l| 
and on the Bethe lattice [2^, extending previous calcu- 
lations Elm of minor loops. 

The equilibrium properties of the RFIM are governed 
by a zero-temperature fixed point, and in finite dimen- 
sions (d < 5 in practice) GS calculations have elucidated 
the properties of the phase diagram. In d > 3 the GS 
displays a ferromagnetic phase transition induced by the 
disorder. As domain wall energy arguments and exact 
mathematical results indicate, in d = 2 there is no phase 
transition but an effective ferromagnetic regime for small 
systems, while in d = 1 the RFIM is trivially paramag- 
netic. It has been suggested that the transition in the GS 
is ruling the transition in the non-equilibrium hysteresis 
loop, also because mean-field calculations give the same 
results in and out of equilibrium psj . Numerical values 
of the exponents are close but not equal, but one must 
consider the difficulties in extrapolating values from the 
finite size scaling [2^, [2^. More recently, the question 
of the universality of the exponents, with respect to the 
shape of the disorder distribution, was discussed in d = 3 
simulations, mean-field theory, and on the Bethe lattice 

Below we report a detailed comparison of the zero 
temperature equilibrium and non-equilibrium properties 
of the RFIM with Gaussian distribution of the random 
fields. We first analyze the problem in d = 1, where exact 
results can be obtained. The average value of the energy 
is computed as a function of the disorder strength for 
the DS and the GS. A direct comparison of the two val- 
ues shows that for weak disorder the differences become 



more substantial, while for strong disorder, where each 
spin basically aligns with the random field, the difference 
tends to vanish. Numerical studies using the same disor- 
der realizations reveal that the main difference between 
the two states comes from the complete reversal of GS 
domains in the DS. This is also visible in the overlap 
between the GS and DS. 

We then study the d = 3 case in which both para- 
magnetic and ferromagnetic behavior exist. The ques- 
tion of whether the transitions appearing in the GS and 
in the hysteresis loop are universal has often been de- 
bated in the literature [28l |2^ | . At the mean-field level 
it is not possible to distinguish the equilibrium and the 
non-equilibrium case and the transition if thus trivially 
the same. In addition, the e expansion for the equilib- 
rium and hysteretic transitions is the same to all orders, 
but one should always consider the possibility of non- 
perturbative corrections to the field theory. Numerical 
simulations in d = 3 indicate that the critical exponents 
and the critical disorder in the two transitions are reason- 
ably close, but the numerical uncertainties do not allow 
for a conclusive statement about their identity. Here we 
directly compare the behavior of the GS and the DS in 
d = 3 close to the disorder induced phase transitions. Wc 
show that while the non universal critical parameter R c 
differs in the two cases, the universal finite-size scaling 
curve for the order parameter can be collapsed on the 
same curve. This suggests some kind of universality in 
the GS and the DS transitions. The numerical simula- 
tions for the GS and DS are done for the same disorder 
realizations for the both cases, for cubic lattices of linear 
sizes L = 10,20,40,80. The results are averaged over 
several realizations of the quenched random fields. In 
both cases, we compute the average magnetization as a 
function of the disorder width. 

A difference in the location of the critical point for 
equilibrium and non-equilibrium behavior of the same 
model may appear rather peculiar and one could be 
tempted to ascribe it to finite size corrections. In order 
to clarify this issue, we have solved exactly the model on 
the Bethe lattice and compared the results for GS and 
DS. While the exponents, as expected, are the same, co- 
inciding with the results of mean-field theory, the critical 
disorder differs in the two cases. Namely the transition 
in the DS occurs at a lower disorder value. Thus there is 
an intermediate parameter region where the GS is ferro- 
magnetic but the DS is paramagnetic. In conclusion, the 
solution on the Bethe lattice corroborates the picture ob- 
tained from simulations in d — 3. From the optimization 
viewpoint, the d = 3 case shows an intermediate phase of 
"bad" correspondence between the GS and DS, exactly 
as in d = 1. This however stops as the Rc is ap- 
proached: naturally if both the states are ferromagnetic 
the optimization of the DS is much easier. To further ex- 
plore the question of universality of the two transitions 
in the GS and in the DS, we have computed the distribu- 
tion of the magnetization at the respective critical point, 
R C DS ^ and r[ GS ^ for different lattice sizes. The distribu- 
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tions can again all be collapsed into the same curve. 

Finally, we consider the question of when is it actually 
possible to reach the exact GS via demagnetization. To 
this end, we consider a reverse field history (RFH) algo- 
rithm that allows in principle to construct a field history 
to get to the GS, if possible at all. Studies of the d = 1 
case illuminate the difficulty of optimizing since it turns 
out that for anything but very strong disorders R the 
probability to reach the GS rapidly decays to zero. 

Our main conclusion is that, in general, demagnetiza- 
tion is not a good technique for reaching states that are 
truly close to the equilibrium, except in cases where the 
outcome is clearly similar from the very beginning (FM 
states and PM states where the disorder is strong) . This 
holds for both the energy of the states and also for the 
spin configurations. A simple formulation is that, since 
the DS is not optimized well in terms of the locations 
of domain walls, it has an excess random field (Zeeman) 
energy. 

The paper is organized as follows: in section II we de- 
fine the model and discuss its numerical treatment. In 
sec. Ill we analyze the one-dimensional case, analytically 
and numerically. Section IV is devoted to the behavior 
around the disorder induced transition in d — 3 and on 
the Bethe lattice. Section V demonstrates the RFH al- 
gorithm, together with numerical studies. Conclusions 
are reported in section VI. An account of some of these 
results was briefly reported in Ref. |34| . 

II. THE RANDOM-FIELD ISING MODEL 

In the RFIM, a spin Sj = ±1 is assigned to each 
site i of a d— dimensional lattice. The spins are cou- 
pled to their nearest-neighbors spins by a ferromagnetic 
interaction of strength J and to the external field H. 
In addition, to each site of the lattice it is associated 
a random field hi taken from a Gaussian probability 
p{h) = exp(-h 2 /2R 2 )/V2nR, with variance R. The 
Hamiltonian thus reads 

n = - Js i s i - J2( H + ft *)**> (*) 

where the first sum is restricted to nearest-neighbors 
pairs. 

In this paper we will consider only the case of zero 
temperature, both in equilibrium and out of equilibrium. 
The T = equilibrium problem amounts to find the mini- 
mum of TL for a given realization of the random- fields (i.e. 
the GS) and then eventually perform the thermodynamic 
limit. This problem has been solved exactly in a num- 
ber of simple cases, namely in d = 1 and on the Bethe 
lattice, for particular disorder distributions and studied 
numerically in generic dimensions. 

The RFIM GS is solvable in a polynomial CPU- 
time, with exact combinatorial algorithms. For the one- 
dimensional case, the solution can be found via a map- 
ping to a "shortest path problem" (3{| which effectively 



places the domain walls in optimal positions, correspond- 
ing to the global minimum of H.. For higher dimensions, 
one starts by noticing that finding the RFIM GS is equiv- 
alent to the min-cut/max-flow problem of combinatorial 
optimization. This can be solved in a variety of ways. 
We use a so-called push-relabel variant of the preflow al- 
gorithm |36j| . Such methods, properly implemented, are 
in general slightly sub-linear in their performance as a 
function of the number of spins in the problem. 

For the out of equilibrium case, we need to specify an 
appropriate dynamics, ruling the evolution of the spins. 
We will consider the dynamics proposed in Ref. [37j and 
used in Refs. 0,0,0 to study the hysteresis loop. At 
each time step the spins align with the local field 

Si = sign( J ^ Sj +hi + H), (2) 

3 

until a metastable state is reached. This dynamics can be 
used to obtain the hysteresis loop. The system is started 
from a state with all the spin down si = — 1 and then H 
is ramped slowly from H — > — oo to H — > oo. The limit 
of dH/dt — > can be conveniently obtained by increasing 
the field precisely of the amount necessary to flip the first 
unstable spin. A single spin flip increases the local field 
of the nearest neighboring spins, generating an avalanche 
of flippings. When the systems finds another metastable 
state, the field is increased again. This dynamics obeys 
return-point memory 0|: if the field is increased adi- 
abatically the magnetization only depends on the state 
in which the field was last reversed. This property has 
been exploited in d = 1 pH |24| and in the Bethe lattice 
[22I |2^ | to obtain exactly the saturation cycle and the 
minor loops. 

The main hysteresis loop selects a series of metastable 
states, which in principle are not particularly close to 
the ground state. To obtain low energy states, we per- 
form a demagnetization procedure: the external field 
is changed through a nested succession H = Hq — > 

Hi -> H 2 -> H n ... -> 0, with H 2n > H 2n+2 > 0, 

H 2n -i < H 2n +\ < and dH = H 2n -H 2n+2 -+ 0. A per- 
fect demagnetization can be performed numerically using 
the prescription discussed above to obtain dH/dt — > 0. 
Such a perfect demagnetization is quite expensive compu- 
tationally and it is convenient to perform an approximate 
demagnetization using dH — 1CU 3 . A comparison of the 
states obtained under approximate and perfect demag- 
netization shows negligible differences. 



III. GROUND STATE AND DEMAGNETIZED 
STATE IN ONE DIMENSION 

A. Exact results: ground state 

The GS energy can be computed exactly in d = 1 using 
transfer matrix methods |9( The free energy of a chain of 
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length N is given by 
Fn = -\ ln(Zjv) = -~ln(Z N - Z N) 



0' 



-~ln(Z+Z N ) 
(3) 

where Zn is the partition function with free boundary 
conditions and Z N are the partition functions with the 
spin at site 1 fixed up (down). These functions satisfy the 
following recursive relation: 



Z^ = e ± ^(Z+_ 1 e^ J + Z^e^ 



j n — c V^JV-1° "T^JV-l ) (4) 
The last step in eq.@ uses the approximation Z^+Z N ~ 



Z N Z N which holds in the large N limit since Z N both 

diverge with the ratio Z N /Z N being finite. From Eq. (0} 
it follows 

Z+Z N = Z+_ 1 Z w _ 1 (2cosh(/3J) + 2cosh(2^jv)) (5) 



where xn = \n(Z~^/Z N ), which gives for the total free 
energy 



Fn = -Fjv-i 



1 

w 



ln(2 cosh(/3 J) + 2 cosh(2/3x w ) ) (6) 



where xn = ^ \n(Zj^/Z N ), so that one can define a free 
energy per site 

/ = --^ln(2cosh/3J + 2cosh(2 ( &r A r)). (7) 

xjv is a stochastic quantity satisfying the equation 

xn = h\ + g(x N -x) (8) 

where g(x) = ± In ((e 2 ^- 7 ) + l) / (e 2 ^) + e W))) 
When i? — > Eq. JSJ has a fixed point solution of 
x oo = g{xoo). It is easy to check that Xqo = is the 
only solution for any J and (3 finite, corresponding to 
the absence of a phase transition. 

When R is non-zero a; at is a random variable with an 
associated distribution Wn(x), where 

Wn{x)cLx = Prob(x < xn < x + dx). (9) 

Wn(x) satisfies the recursive functional equation 

W N+1 (x) = J^dhP^x 

J^dxiWNfaMx-h-H-gfx!)) (10) 

so that in the thermodynamic limit is given by the 
fixed point equation 



vMx) 



dxiWoo(xi)P(x - h - H - g(xi)) . (11) 



Once Wqo is known, any thermodynamic quantity can be 
computed. In particular, the free energy per spin, which 
is given by 

1 f°° 

(/) = -- / dxWootx) (cosh(2/3) + cosh2/3x) . (12) 

PJ-oo 



The magnetization at a site of an infinite lattice, is 
given by 

to \ - zT-z-t _ 
\i>0/ — zT+Zl ~ 

^^ = tanh(Il n( ZVZi)), (13) 

where Z'^ are respectively the partition functions with 
the spin at fixed up (down). These are given by 



Z^: 



Z~)(e 



±/3 J 7+ 



(14) 

where Z r l are the partition functions for the semi-infinite 
right(left) lattice, with the spin at site 1 (—1) fixed 
up(down). This gives 



(s ) = tanh(/3(/i + g(x r ) + g(xi)))- 



(15) 



Finally, The magnetization for the infinite lattice is ob- 
tained averaging over the quenched variables x r ^i : 

m = j 00 ^ dhP(h) dx r W N {x r ) 
d Xl W N (xi) tanh [fi{h + g{x r ) + g( Xl ))) . (16) 



B. Exact results: Demagnetized state 

In d = 1 the magnetization and the energy per spin 
as a function of the external field can be derived explic- 
itly through a probabilistic reasoning. We show how to 
get these results on the saturation loop, focusing on the 
lower branch. (The results on the upper branch can be 
obtained by symmetry considerations.) Similar but much 
more involved reasoning can be repeated for any minor 
loop. 

The central quantity to consider, in order to solve for 
the magnetization as a function of the external field H 
on the hysteresis loop, is the conditional probability for a 
spin to be up, conditioned to one of its nearest neighbors 
being down. To calculate this quantity, one can reason 
as follows: fix the spin at site i — 1 down. Define p m (H) 
as the probability for a spin to be up, given that exactly 
m (m = 0, 1, 2) of its neighbors are up: 



p m (H) = P{hf > 0) = J 



dhp(h) , (17) 



- 2m)./ - H 



where z is the coordination number (z = 2 in d = 1). Fix 
for a moment the spin at site i down as well and look at 
the spin at site i + X. It will be up with probability Uq 
and down with probability 1 — Uq. The spin at site i will 
flip up with probability pi when the spin at i + 1 is up 
and po when it is down. Ultimately, the spin at i will 
be up (conditioned to the spin at i— 1 being down) with 
probability Uq = Uopi + (1 — f/o)po- It follows 



Uq 



Po 



1 - Pi + Po 



(18) 
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Once Uq is known, a similar reasoning leads to the (un- 
conditioned) probability p(H) for a spin to be up: Fix 
the spin at site i down. The spin at site i — 1 will be up 
with probability Uq and down with probability 1 — Uq. 
The same holds for the spin at site Thus 

P (H) = U 2 p 2 + 2U (1 - U )P! + (1 - U ) 2 P(h (19) 

from which the magnetization is obtained as m(H) = 
2p(H) - 1. 
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FIG. 1: The energy of the GS is compared with the one of 
DS. The values are computed exactly in d = 1 as a function 
of the disorder width R. 
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FIG. 2: The energy difference between the GS and the DS 
computed exactly in d — 1 as a function of the disorder width 
R. 

The energy per spin on the saturation loop is obtained 
as follows. Due to translational invariance: 



E 



N 



-J(si8i + t) - His,,} - (hiSi). 



(20) 



To calculate the spin-spin correlation (s,-Sj_)_i) we intro- 
duce the probabilities <I> ++ , $ + ~, $~ + , for adja- 
cent spins to be respectively up-up, up-down, down- 
up, and down-down. These quantities are not indepen- 
dent, since they have to satisfy the obvious identities: 
$+-=$-+, $+++$+- = p(H), and*— +<P + - = l-p(H). 
Thus it is sufficient to calculate one of them, for ex- 
ample This is done by separating the four con- 
tributions from the possible boundary conditions deter- 
mined by the values of the spins at sites i — 1 and i + 2: 
When they are both down, the probability for the cou- 
ple of spins at sites i and i + 1 to be both down is 
Uq (1— pi(H)) 2 , when one is up and the other one is down 
it is 2f7 (l - E7o)(l -pi(H))(l -p (H)), and when both 
of them are up it is (1 — Uq) 2 (1 — pq(H)) 2 . Adding up the 
four contributions one gets = (1 — Uq) 2 . This fixes 
the other probabilities to be $ + ~ = <f>~ + = 2p-l+(l-C/ ) 2 , 
and <f> ++ = 1— p — (1 — Uq] 2 . Thus, the spin-spin corre- 
lation is 



{ Sl s l+1 ) = $+++ 



2$ 4 



A{p-{\~Uq) 2 ) 



(21) 

The average value (hiSi) can be obtained by averaging 
over the held h the product of h times the average value 
of the spin Sj over the local fields other then hi , once the 
field at i is fixed at the value h : 



(hiSi) 



dh p(h )h (Si\h ). 



(22) 



The conditional average (si\h ) is given by 2p(H\h') — l 
where p(H\h ) is the conditional probability for a spin 
to be up at an external field H , given that its local ran- 
dom field is fixed at the value h . From Eq. (|19fl this is 
trivially given by 



p(H\h) = U%6{ti + H + 2 J) 

+ 2U (1 - U )6(h + H) 
+ (1 - U ) 2 e(h + H - 2 J) 

which finally gives 

<> + oo 



(23) 



{hiSi) = 2U 2 



dh p(h )h 



+ Wq{\-Uq) 

+ 2{1-U ) 2 



dh p(h )h 

dhp{h)h~h'. (24) 



In particular, for a Gaussian distribution with h' — and 
variance R the integrals can be performed analytically 
and the result is 



2 r 

Re 2r 

7T 



2U^eW cosh (2JH/R 2 



+ e 2J(J -^ ) (l-2U 2 ) + 2UQ(l-U )\ .(25) 
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The energy per site on the lower branch of the saturation 
loop is in general given by 

E(H) = -47 (p(H) - (1 - U ) 2 ) + 3 J - H(2 P (H) - 1) 

/+oo 
dh p{h)h 
H — 2J 

/+oo 
dh p{h)h 
H 



+ 2(1 - U )' z I dh p(h )h - K . 

-H + 2J 



(26) 



Similar but much more involved reasonings can be re- 
peated for any minor loop - eventually for a series of 
nested loops leading to the demagnetized state - pro- 
viding a series of recursive equations for the magneti- 
zation, the spin-spin, and the spin-field correlations, 
which are the quantities needed to compute the en- 
ergy. If the external field is changed through a nested 

succession H = H — > Hi — > H 2 — ► H n ... — > 0, 

with H 2n > H 2n+2 > 0, H 2n -i < H 2n+ i < and 
dH = H 2n — H 2n+2 — > 0, the spin-spin correlations are 
given recursively by 



4Ui n ( P2 (H 2n )-p 2 (H 2n ^)) 

Po(H 2n -i)) 



(27) 



where Uj. and Dk are respectively the probabilities for 
a spin to be up (down) conditioned to one of its neigh- 
bors being down, and satisfy in turn a set of recursive 
equations. Similar equations hold for magnetization and 
spin-field correlation, leading to a complicated recursive 
formula for the energy. The results of such calculations 
are shown in Figs. (QEJ), where the energy of the demag- 
netized state is compared with the energy of the ground 
state evaluated in the previous section. 



C. Simulations: how optimized is the 
demagnetized state? 

In one dimension the comparison of the DS and the 
GS is the easiest since the domain walls are just point- 
like. For the GS we know that it is optimized such that 
all the large enough local random field fluctuations nu- 
cleate domains of the same sign. The rest of the random 
landscape is split up into regions that align themselves 
with such fluctuations depending on the sign of the ran- 
dom field excess, XiGrcgion hi- A s a resu lt the Zeeman 
energy of domains is linear in domain size, Ez ~ Id, and 
the asymptotic mean domain length follows the Imry-Ma 
prediction (Igs) ~ 1/R 2 - Moreover since the random 
landscape has a finite correlation length the domain size 
distribution is exponential j3^. 

Any qualitative differences in the DS will follow from 
three separate mechanisms: 1) shifts of domain walls, 
2) creation of domains inside intact GS domains and 3) 



destruction of GS domains (Fig. UJ . From the point of 
view of "optimization" the first one is of trivial concern, 
since it would have little effect e.g. on the scaling of 
Ez,ds- The second one is more detrimental if the energy 
difference to the GS is considered. In addition to the 
cost of the two domain walls it subtracts a contribution 
from the Zeeman energy of the domain that persists and 
surrounds (in the GS) the one that is not created in the 
DS. The third one would make the largest change to the 
total energy, since for Id 3> 1 the energy of a domain 
consists mostly of its Zeeman energy. 

+ + + - - -- -- + + + groundstate 

, 1 

+ --+++ DW shift 

i_ J 

+ + + - -'++!- - + + + nucleated 
i ■ 

. . . ~ ~ ; . droplet 
+ + + + + + + ;+ + + destroyed 

FIG. 3: An illustration of the possible mechanisms for the 
deviations between GS and DS. 
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FIG. 4: The average change in the spin-spin overlap between 
the GS and the DS (Aq) and the contribution to that from 
completely "destroyed" GS domains (Ag destr .), as a function 
of R. 

Numerical studies of the DS domain structure indicate 
that with decreasing R the average domain size increases 
faster than in the GS, while the size distribution P(ld) 
remains exponential. This is accompanied by a reduction 
in the overlap q = ((ugso'ds) + l)/2 between these two 
states. For R large the overlap is close to unity; strong 
local fields hi align the spins in the same way regardless 
of the mechanism by which the spin state is created. For 
R small the local field is no longer strongly correlated 
with the orientation of the spin, and thus whether the 
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0.6, 0.8, 0.9, 1.0. The overlap decreases with Idomain.DS ■ 
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FIG. 6: The Zeeman energy of DS domains as a function of 
the size, ldomain,DS- The black circles mark the average DS 
domain size for a given R. The two lines above and below 
the data indicate optimal, linear (GS-like) scaling and the 
Imry-Ma-like /^-scaling, respectively. 



GS and DS are locally aligned depends on how optimized 
the latter is. 

The fundamental mechanism for the deviations be- 
tween the states seems for R small to be the "destruction" 
of GS domains (see Fig. [3] again) . This is demonstrated 
in Fig. 2] by depicting the change Aq in the overlap that 
comes solely from missing GS domains. The conclusion 
from this dominance is that the demagnetized states typi- 
cally miss regions in which the integrated field fluctuation 
is large which as such leads in the GS to the formation of 
GS domain. Therefore the overlap should get smaller the 
larger the scale-length on which one compares the DS and 
GS is, and this is confirmed by Fig. [S] which shows the 



overlap between a DS domain and the GS as a function 
of the length of the DS domain. 

The importance of such destroyed domains can also 
be seen in the total contribution to the energy difference 
between the DS and GS. For R small this is again domi- 
nated by the missing GS domains. In general the differ- 
ence between the energies of the GS and DS derives from 
the combination of domain walls and Zeeman energy. 
Fig. [B] shows that for Id small the DS domains do not 
have much Zeeman energy. This changes if Id is larger, 
in which regime the scaling approaches the Imry-Ma -like 
scaling (7°' 5 ). The implication is that the field energy of 
large domains in the DS self-averages, and comes from a 
sum of random contributions (ie. the domains contain re- 
gions where the actual random field sum is opposite to the 
spin orientation, such as the missing GS domains). The 
cross-over between the small /^-behavior and the asymp- 
totic scaling is located close to (ld)DS- 

IV. AROUND THE DISORDER INDUCED 
TRANSITION 

A. Simulations in d — 3 

The RFIM displays a disorder induced phase transition 
both in the GS and in the hysteresis loop, which can also 
be observed analyzing the DS 0, |2^, H^J . If the GS and 
the DS are always paramagnetic, the transition is absent 
and thus we perform numerical simulations in d = 3. Our 
aim is to characterize the difference between DS and GS 
around the disorder induced transition. 

In d = 3 for low disorder, the GS is ferromagnetic, 
while for higher disorder it becomes paramagnetic. For 
Gaussian disorder, the transition point has been located 
numerically at i?c GS ^ — 2.28. It is possible to define the 
usual set of critical exponents characterizing the phase 
transition and compute the values by exact GS calcula- 
tions. For instance, the magnetization M = (\m\), with 
m = J^. Si/N, scales close to the transition point as 

M = Ar?, (28) 

where r = (R — R c )/R c is the reduced order parame- 
ter and A is a non-universal constant. The correlation 
length defines another exponent £ = (Br)^ 1 ' -where B 
is another non-universal constant- which rules the finite 
size scaling of the model 

M = AL-M v f (BL^^R - R c )/R^j . (29) 

Simulations yield ~ 1.17 and /3 (GS) = 0.02. 

A disorder induced transition is also found in the hys- 
teresis loop. At low disorder the loop shows a macro- 
scopic jump, which disappears at a critical value for the 
disorder. This transition reflects itself in the DS, which 
is ferromagnetic when the main loop has a jump and is 
paramagnetic otherwise. The transition point has been 
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i 2.16 and the 
In particular, 



obtained numerically in d — 3 as R^P S ^ 
critical exponents have been measured. 
Ref. [3£| reports data collapses with /3ms) — 0-04 and 
v (ds) — 1-41. While there is strong evidence that the 
exponents measured in the DS should be equal to those 
measured on the main loop, the relation with the equi- 
librium transition is not clear. 

We notice first that numerical simulations reported 
in the literature indicate that the transition appears at 
slightly different locations in the GS and in the DS. Hart- 
mann and Nowak report R^ s = 2.29 ± 0.04 for the GS 
with system size up to L = 80, Hartmann and Young 
refine this value to Rc GS ^ = 2.28 ± 0.01 with sizes up to 

L = 96, which is also confirmed by Middleton and Fisher- 
ies'') 

which estimate R c = 2.27 ± 0.04. For the hysteresis 
loop the best estimate is R c — 2.16 ± 0.03, with sys- 
tem sizes up to L = 320 and a similar value for the DS 
0,|3^|. Thus, unless strong finite size effects take place, 
one is tempted to conclude that the two transitions take 
place at two different values of R c . 
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FIG. 8: Numerical results in d = 3: The magnetization can be 
collapsed using R c = 2.28 (GS) and R c = 2.16 (DS), l/v = 
0.73 and /3 = 0.03. The scaling curve is the same for DS 
and GS indicating universal behavior. The values for the 
ratios of the non- universal constants are Ads /Acs = 1 and 
B D s/Bgs = 0.68. 



FIG. 7: The magnetization of the GS and the DS in d - 
different system sizes L and disorder R. 



3 for 




Here we analyze the problem again by numerical simu- 
lations, computing the GS and the DS numerically, using 
the same disorder realizations for the two cases. Sim- 
ulations are performed for cubic lattices of linear sizes 
L = 10,20,40,60,80 and the results are averaged over 
several realizations of the random fields. The GS is found 
exactly using a min- cut /max- flow algorithm, while de- 
magnetization is performed approximately with the algo- 
rithm discussed in Ref. [2l| with dH = 10 -3 (see section 
II). In both cases, we compute the average magnetiza- 
tion as a function of the disorder width (see Fig. 01. In 
Fig. [S] we collapse the two sets of data into a single curve, 
using two different values for R c (i.e. R c — 2.28 and 
r[ ds ^ = 2.16) but the same values for the exponents (i.e. 
1/v = 0.73 and (3 — 0.03). The best value for the ratio of 
the non- universal constant is found to be Ads/^-gs — 1 



FIG. 9: The overlap between the GS and the DS in d ■ 
different system sizes. 



3 for 



and B DS /B GS = 0.68 ± 0.02. The fact that the scal- 
ing function is the same for the two cases is a strong 
indication for universality, going beyond the simple nu- 
merical similarity of the exponents. There is always the 
possibility that in the limit L — > 00 i£c G = Ri DS ^ . At 
the present stage this hypothesis is not supported by the 
data, since we were not able to collapse all the data into 
a single curve using the same R c . 

Next, we compare the statistical properties of the GS 
and the DS around the transitions. In Fig. we report 
the value of the overlap as a function of R for different 
system sizes. When the disorder is decreased from the 
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FIG. 10: The difference in magnetization between the GS and 
the DS in d = 3 for different system sizes. 



paramagnetic region, the overlap decreases as for d = 1. 
However for low disorder the overlap rapidly increases 
and reaches 1 in the ferromagnetic state. The minimum 
of the overlap is located in the parameter region cor- 
responding to the transitions (i.e. R ~ 2.2 — 2.3). A 
decrease in the overlap around the transition can be ex- 
pected, since for Rc S ' < R < Ri GS ^ the GS is ferro- 
magnetic (M > 0) and the DS is paramagnetic (M = 0) 
as it is also apparent plotting the difference in the mag- 
netization (see Fig. lTfj|) . 

In summary, three dimensional simulations indicate 
that the transitions in the GS and DS are universal, but 
the critical parameter seems to differ. Consequently the 
GS and DS differ mostly around the transition, while the 
difference is smaller in the paramagnetic and ferromag- 
netic phases. 



B. The Bethe lattice 



where 



1 



Fn-i{j)-7^ In 2 (cosh(/3J) + cosh(2/3:E n (j))) 

(31) 



x n {i) = —Mz+(i)/z-(i)), 



(32) 



so that the contribution at the free energy from site i is 
f(i) = -— ln(2cosh0J + 2cosh(2/?ir„(i))). (33) 
x n (i) is a stochastic quantity satisfying the equation 

X n (i)=hi+ 2J 9( x n-l(j)) ( 34 ) 

jei(i) 

When R — > Eq. I|34|) has a fixed point solution of 
= (z — l)g(x oa ). — is a solution for any J 

and p. For /3 < (3 C = | In there are also two stable 
solutions ixoo 5^ corresponding to the appearance of a 
ferromagnetic phase. 
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The RFIM can be solved exactly in the Bethe lattice, 
displaying a disorder induced transition in the GS and in 
the DS [22|. It is thus an interesting case to compare the 
two states around the respective transition directly in the 
thermodynamic limit. We consider here a Bethe lattice 
with coordination z and obtain the GS generalizing the 
d = 1 case as in Ref. |38j. In this case N refers to 
the generation of the lattice, and are the partition 
functions of a branch of generation n with a fixed up 
(down) spin at the central site. The recursion relation 
for the Z^ is 

Z±(i)=e ± ^ J] (Z+_ 1 U)^ J + Z-_ 1 (j)e^ J ) (30) 

where for any given site i the sum over j runs over the 
z — 1 nearest neighbors of i away from the center of the 



FIG. 11: The magnetization of the GS and the DS computed 
exactly on the Bethe lattice with z — 4 in the thermodynamic 
limit, showing the ordering of the critical point (see inset). 
When the data are plotted against the reduced parameter 
(R c — R)/R c the curves superimpose. The result implies that 
for the Bethe lattice Acs = Ads- 

To perform quenched averages one has to solve for the 
probability distribution of W n (x n ), where W n {x)dx = 
Prob(x < x n < x + dx), which satisfies the recursive 
functional equation 

/oc />oo 
dhP(h) / d Xl W n {xx) ■ ■ ■ 
-oo J — oo 

/CO f ~3 

dx z ^W n (x z ^)6(x -h~H -E<?(zfc)) .(35) 
-oo 
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so that in the thermodynamic limit Woo is given by the 
fixed point equation 



Woo 0*0 = / dxiWoo(xi) ■ ■ ■ 

J —OO 

dx z .W oix z - 1 )P{x- h- H -^ 5 (x fe )).(36) 
-°° fc=i, 

Once Woo is known, any thermodynamic quantity can 
be computed. In particular, the free energy per spin is 
given again by l|12|l and the magnetization at the central 
site of an infinite lattice, is given by Eq. I|13l) where 
are respectively the partition function with the spin at 
fixed up (down). These are given by 



z n = e ±0h o 



J] (e^t + e^-) (37) 



k=l, 



and Zj^ for k 



1, • • • z are the partition functions of the 
z branches attached to the central site 0, with the bound- 
ary spin fixed up(down). This gives for the magnetization 
at the central site (s ) 



(s ) = tanh(/3(/i + ^ g(x k ))) 



(38) 



fc=i, 



The magnetization for the infinite lattice can then be 
obtained averaging over the quenched variables x r j: 



M I dhP(h) / dxiW N (x{)--- 

) J — OO 

dx z W N (x z )tanh(/3(h+ ( 39 ) 



fe=i, 



For a Gaussian random-field distribution the fixed point 
equation can not be solved explicitly and we thus re- 
sort to a numerical integration. We obtain Woo{x) for 
z = 4, and for different values of R, and compute the 
magnetization using Eq. Ij39(l . In Fig. 1111 we compare 
the magnetization of the GS with the one of the rem- 
nant magnetization in the DS, computed in Ref. j^. As 
observed in the simulations in d = 3, the transition oc- 
curs at two different locations (see the inset of Fig. Ill|l . 



for z = 4 R { C DS) = 1.781258... [2S 



and R, 



(GS) 



1.8375, 



with the mean- field exponent {(3 — 1/2). When plotted 
against (R — R c )/R c the two curves superimpose close 
to the critical point. This indicates that, though not re- 
quired by universality, in the Bethe lattice Aqs — Ajjs, 
as also found in d = 3. 

To investigate possible finite size scaling we have per- 
formed numerical simulations in the Bethe lattice, fol- 
lowing the method of Ref. |25j . Collapsing the order pa- 
rameter curve as in d — 3, using a scaling form similar 
to Eq. J2|J), does not appear to be possible in the Bethe 
lattice, because the scaling region is very narrow. Thus 
to test finite size scaling, we have computed the distri- 
bution of the magnetization m at the respective critical 



point, i?c DS ' ) and Ri ^ for different lattice sizes N. The 
distributions can all be collapsed into the same curve (see 
Fig.Q2J}, using the form P(\m\) = f(\m\/M)/M. 
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FIG. 12: The distributions of the magnetization in the DS 
and the GS at their respective critical points on the Bethe 
lattice, obtained numerically for different lattice sizes iV, can 
be all collapsed together. 



V. REACHING THE GROUND STATE BY 
NON-EQUILIBRIUM DYNAMICS 

After having shown that the GS and the DS corre- 
spond to different microscopic configurations, we investi- 
gate now if the GS spin configuration may be reached by 
a field history other then the ac-demagnetization. The 
answer to this question requires a clarification on the 
relation existing between locally stable states (given by 
the solutions of Eq.Q), and the spin configurations vis- 
ited along the non-equilibrium dynamics induced by the 
varying field. In fact, not all stable configurations may 
be reached by a field history from saturation. The prob- 
lem has been treated in |40j where it has been shown 
that, given a spin configuration obtained by a field his- 
tory, supposed unknown, the sequence of reversal fields 
that applied to saturation gives back the original state 
can be recovered. For spin systems this inverse function 
is given by an algorithm which is able to construct the 
reverse field history (RFH) ^(|- This method is applied 
then to investigate if a given spin configuration may be 
reached by field history from saturation: if a field his- 
tory leading to the state exists the algorithm produce a 
sequence of reversal fields; if no field history exists the 
algorithm enters a recursive loop. The investigation of 
the properties of unreachable states has been recently 
performed and leads to a classification of stable states 
on oriented graphs ^lj- The study is performed here for 
the GS spin configuration that, for the RFIM at finite 
size and for a given disorder realization, can be indepen- 
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dently derived by exact combinatorial algorithms (as the 
max- flow min-cut). 



A. RFH Algorithm 

Consider the final spin configuration s (the set of N 
Ising spins) resulting after the application of a field his- 
tory ending at H = and consisting in a sequence of 
reversal fields {H} = {Hi, . . . , H n } from the saturated 
state and let us define the function s = f({H}). The set 
of all states obtained this way is defined as the hysteresis 
states (H-states). Due to adiabatic dynamical response 
and return point memory, this state s will contain the 
memory of a subset of the reversal fields. In fact not all 
the reversal fields determine the final state s. For ex- 
ample, in terms of average magnetization, the reversal 
fields which give rise to closed minor loops do not influ- 
ence the final state, i.e. their memory is erased, while 
the memory of the set of reversal fields {Hs} which are 
not erased is contained in the final state. The inverse 
function {Hs} = g(s) allows, starting from a spin config- 
uration s belonging to the H-states ensemble, to obtain 
the set of reversal fields {Hs} which have been actually 
stored in the state and that - if applied as a field history 
- will reproduce the original state, i.e. s = f(g(s)). We 
define this set of reversal fields {Hs} as minimal field 
history. 

The RFH algorithm takes as input a configuration s at 
H = and gives as output - when it exists - the reversal 
field history from saturation to the state s. The formu- 
lation of the algorithm is based on the order-preserving 
character of the dynamics [l^, and is therefore, appli- 
cable to a wide range of models beyond the RFIM. An 
interesting result of the RFH algorithm is obtained when 
it is applied to a state s not belonging to the H-states 
(i.e. where no field history exists). The iterated search 
for the reversal field sequence enters an iteration and, in 
this case, it can be shown that no field history leading to 
the state exists. 



As a first finding the GS does not result to be sys- 
tematically field reachable and the fraction depends on 
the disorder. One may conclude that the fact that the 
GS is sometimes reachable is a pure effect of the finite 
system size. However, also for the random states the 
fraction of found states ijiND sensibly changes with R, 
but following a different curve. If there was no correla- 
tion between GS and the H-states the two curves would 
be coincident. The dependence of fRND on R reflects the 
fact that the number of H-states depends on the disorder 
value and the system size |4^ . and only at large disorder, 
where the number of locally stable states decreases, the 
ratio between H-states and stable states is significantly 
greater then zero. 




FIG. 13: Fraction of reachable states (averages over 100 real- 
izations of disorder) diamonds: fraction Jrnd ; circles: frac- 
tion fas 



VI. CONCLUSIONS 



B. Simulation results in Id 

The RFH algorithm was applied to explore the possi- 
bility to reach the GS by non equilibrium dynamics by 
the numerical study of the RFIM in one dimension with 
periodic boundary conditions. We performed our investi- 
gations on systems having N — 5000 spins, averaging the 
results for 100 different realizations of the same disorder 
R. The GS was obtained by the max-flow min-cut pro- 
cedure for each realization and the RFH algorithm was 
applied. At each disorder value R, the fraction fas of 
the realizations in which the GS resulted to be reachable 
was computed. For comparison the same procedure was 
applied starting from locally stable states generated by 
random sampling the set of local minima. The results 
are shown in Fig lloi 



For disordered systems like the random field Ising 
model one would be interested in both universality in sta- 
tistical properties and in the question how to "optimize" 
in the case of a sample with a given distribution of the 
impurities. In this paper we have studied this problem 
in detail, by comparing the demagnetized and ground 
states. Our main findings are the following: First, the 
character of the GS is such that it is globally optimized, 
and the demagnetization procedure does not perform well 
unless the optimization problem is rather trivial. This is 
slightly surprising since the conclusion holds in particu- 
lar if the RFIM GS is paramagnetic. Then the DS does 
not manage to find the right spin configuration, so that 
as seen in the d = 1 case many of the domains of the GS 
do not appear in the DS. 

Second, in d = 3 (and with the aid of the Bethe 
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lattice solution), it can be demonstrated that the ex- 
istence of a phase transition for both the DS and GS 
makes the "phase diagram" of optimization to show a 
regime where the outcome is less optimal: in the para- 
magnetic phase of the DS, where the GS is already fer- 
romagnetic since the critical thresholds are ordered such 
that R { C GS) > R ( C DS) = 1.84. In this regime DS and 
GS are expected to differ strongly in the thermodynamic 
limit. We also provide numerical evidence that the d = 3 
transition appears to have the same critical exponents in 
both the GS and DS This can be considered both 

surprising - there being no exact field theoretical way 
of treating the d = 3 phase transition - and expected, 
since the functional renormalization calculations in spite 
of their shortcomings indicate that the actions are the 
same jlij . It seems intriguing that such universality is 
met exactly in the limit where the "optimized" character 
of the DS changes. 

The results indicate that for the particular system at 
hand, where the disorder couples directly to the expected 
magnetization, "local" optimization methods have diffi- 
culties. Of course, as in "hysteretic optimization" , one 
can perturb or "shake" the state obtained from the DS 
procedure to try to still lower the energy. These attempts 
are of course usually just heuristic. In the case of the 



RFIM, the joint approach of optimizing by the DS and 
computing the GS exactly allows to understand better 
similarities and differences between equilibrium and low 
energy non-equilibrium states. 

In addition to the ferromagnetic RFIM model, one can 
consider other systems where two disorder induced phase 
transitions exist. Numerical simulations and analytical 
results have shown that a disorder induced transition in 
the hysteresis loop can be observed in the random bond 
Ising model in the random field O(N) model in 
the random anisotropy model 14(1 147| and in the random 
Blume- Emery-Griffith model |44|. All these systems dis- 
play as well a transition in equilibrium and it would be 
interesting to compare their DS and GS. 

Interfaces in quenched disorder would provide another 
interesting example, since the roughness exponent typi- 
cally differs in and out of equilibrium (i.e. at the depin- 
ning threshold) 4J. It would be interesting to measure 
the roughness of an interface after a demagnetization cy- 
cle (i.e. after the field driving the interface is cycled with 
decreasing amplitude), and compare its properties with 
those of the ground state interface. Finally, there is the 
issue of energetics of excitations in the respective ensem- 
bles: the universality of exponents and scaling functions 
would seem to imply that these also scale similarly. 
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We compare the ground state of the random-field Ising model with Gaussian distributed random 
fields, with its non-equilibrium hysteretic counterpart, the demagnetized state. This is a low energy 
state obtained by a sequence of slow magnetic field oscillations with decreasing amplitude. The main 
concern is how optimized the demagnetized state is with respect to the best-possible ground state. 
Exact results for the energy in d = 1 show that in a paramagnet, with finite spin-spin correlations, 
there is a significant difference in the energies if the disorder is not so strong that the states are 
trivially almost alike. We use numerical simulations to better characterize the difference between 
the ground state and the demagnetized state. For d > 3 the random-field Ising model displays a 
disorder induced phase transition between a paramagnetic and a ferromagnetic state. The locations 
of the critical points Ri DS \ Ri° S ^ differ for the demagnetized state and ground state. Consequently, 
it is in this regime that the optimization of the demagnetized stat is the worst whereas both deep 
in the paramagnetic regime and in the ferromagnetic one the states resemble each other to a great 
extent. We argue based on the numerics that in d = 3 the scaling at the transition is the same in 
the demagnetized and ground states. This claim is corroborated by the exact solution of the model 
on the Bethe lattice, where the R c 's are also different. 



PACS numbers: 

I. INTRODUCTION 

The relation between equilibrium and non-equilibrium 
states is a central problem in the physics of disordered 
systems. Disorder induces a multitude of metastable 
states in which the system can easily be trapped. The 
dynamics is usually very slow, or glassy, and on obser- 
vational timescales the system is basically always out of 
equilibrium. On the other hand, from the theoretical 
point of view it is easier to consider equilibrium proper- 
ties, since in this case it is possible to use all the ma- 
chinery of statistical physics to tackle the problem. The 
question is whether the equilibrium properties of disor- 
dered systems provide a faithful representation of the 
non-equilibrium states in which the system is likely to 
be found in practice. This dichotomy is at the core of 
many unsolved issues in the field of disordered system. 
Typical quantities that one could compare are the energy, 
the geometric characterization of the state (as domains 
in magnets), and the energy cost of excitations. 

A simplification of the problem is obtained considering 
only athermal processes, in which the temperature of the 
system plays no role and can be ignored. The equilib- 
rium state is in this case just the ground state (GS), the 
state of minimal energy [1]. A zero temperature, non - 
equilibrium dynamics is purely relaxational: the system 
falls simply in the closest metastable state. A convenient 
way to allow the system to explore the various metastable 
states is by applying an external magnetic field. Differ- 
ent field histories typically result in hysteresis and lead 



to different metastable configurations [2]. 

The demagnetization process consists in applying a 
slowly varying AC field with decreasing amplitude, and 
provides a simple way to access low energy states [2]. 
It has been studied for more than a century, but until 
recently the question how close the demagnetized state 
(DS) is to the true GS was not addressed. This is the 
concern of our work, the problem of how such an opti- 
mization process works in the case of a random magnet. 
Recently, Pazmandy et al. have proposed the demag- 
netization process as the basis for a new optimization 
algorithm for disordered systems [3]. The idea behind 
such "hysteretic optimization", is that demagnetization 
leads to a low energy state, sufficiently close to the GS, 
which can then be reached by applying other methods 
using the DS as an input. The method was tested for 
different models like spin glasses and NP-hard problems. 

Here wc will concentrate on the random-field Ising 
model (RFIM), that, while retaining some complex fea- 
tures characteristic of disordered systems, still allows for 
a theoretical analysis [4]. In the RFIM, due to the ab- 
sence of frustration, the equilibrium state is relatively 
simple, however, the non-equilibrium dynamics is far 
from trivial. Due to the coupling of the local disorder 
to the order parameter, even the GS presents a variety 
of phenomena, which can be studied numerically [5-8]. 
In fact the GS of the RFIM can be found in a poly- 
nomial CPU-time, with exact combinatorial algorithms 
[1] and solved exactly in d = 1 and on the Bethe lat- 
tice [9, 10]. The equilibrium critical exponents for ran- 



dom field magnets have been measured experimentally in 
Feo.93Zno.07F2 [11, 12]. 

The hysteretic properties of non equilibrium RFIM 
have been widely studied in the recent literature. The 
hysteresis loops display a disorder induced phase transi- 
tion: for low disorder the loop has a macroscopic jump 
at the coercive field, while at high disorder the loop is 
smooth, at least on the macroscopic scale [13-15]. At 
smaller scale the magnetization curve is highly discon- 
tinuous, showing Barkhausen-type bursts, in correspon- 
dence to jumps between different metastable states [16]. 
A disorder induced non-equilibrium phase transition in 
the hysteresis loop has been studied experimentally in 
C0-C0O films [17] and Cu-Al-Mn alloys [18]. 

Extensive numerical simulations have been used to 
characterize disorder induced transitions in the non- 
equilibrium RFIM and critical exponents have been es- 
timated in several dimensions [15, 19, 20]. The model 
has been studied by the renormalization group and the 
exponents have been computed in a e = 6 — d expan- 
sion [14]. In addition the hysteresis loop has been com- 
puted exactly in d = 1 and on the Bcthe lattice, where 
the disorder induced transition is present for sufficiently 
large coordination number. While in d = 1 there is def- 
initely no transition, the situation in d = 2 is less clear. 
Recently the problem of minor loops has been tackled 
analytically and numerically. In particular, the demag- 
netization curve has been computed exactly in d = 1 [21] 
and on the Bethe lattice [22], extending previous calcu- 
lations [23-26] of minor loops. 

The equilibrium properties of the RFIM are governed 
by a zero-temperature fixed point, and in finite dimen- 
sions (d < 5 in practice) GS calculations have elucidated 
the properties of the phase diagram. In d > 3 the GS 
displays a ferromagnetic phase transition induced by the 
disorder. As domain wall energy arguments and exact 
mathematical results indicate, in d = 2 there is no phase 
transition but an effective ferromagnetic regime for small 
systems, while in d = 1 the RFIM is trivially paramag- 
netic. It has been suggested that the transition in the GS 
is ruling the transition in the non-equilibrium hysteresis 
loop, also because mean-field calculations give the same 
results in and out of equilibrium [28]. Numerical values 
of the exponents are close but not equal, but one must 
consider the difficulties in extrapolating values from the 
finite size scaling [28, 29]. More recently, the question 
of the universality of the exponents, with respect to the 
shape of the disorder distribution, was discussed in d = 3 
simulations, mean-field theory, and on the Bethe lattice 
[30-32]. 

Below we report a detailed comparison of the zero 
temperature equilibrium and non-equilibrium properties 
of the RFIM with Gaussian distribution of the random 
fields. We first analyze the problem in d — 1, where exact 
results can be obtained. The average value of the energy 
is computed as a function of the disorder strength for 
the DS and the GS. A direct comparison of the two val- 
ues shows that for weak disorder the differences become 



more substantial, while for strong disorder, where each 
spin basically aligns with the random field, the difference 
tends to vanish. Numerical studies using the same disor- 
der realizations reveal that the main difference between 
the two states comes from the complete reversal of GS 
domains in the DS. This is also visible in the overlap 
between the GS and DS. 

We then study the d — 3 case in which both para- 
magnetic and ferromagnetic behavior exist. The ques- 
tion of whether the transitions appearing in the GS and 
in the hysteresis loop are universal has often been de- 
bated in the literature [28, 29]. At the mean-field level 
it is not possible to distinguish the equilibrium and the 
non-equilibrium case and the transition if thus trivially 
the same. In addition, the e expansion for the equilib- 
rium and hysteretic transitions is the same to all orders, 
but one should always consider the possibility of non- 
perturbative corrections to the field theory. Numerical 
simulations in d = 3 indicate that the critical exponents 
and the critical disorder in the two transitions are reason- 
ably close, but the numerical uncertainties do not allow 
for a conclusive statement about their identity. Here we 
directly compare the behavior of the GS and the DS in 
d — 3 close to the disorder induced phase transitions. We 
show that while the non universal critical parameter R c 
differs in the two cases, the universal finite-size scaling 
curve for the order parameter can be collapsed on the 
same curve. This suggests some kind of universality in 
the GS and the DS transitions. The numerical simula- 
tions for the GS and DS are done for the same disorder 
realizations for the both cases, for cubic lattices of linear 
sizes L — 10,20,40,80. The results are averaged over 
several realizations of the quenched random fields. In 
both cases, we compute the average magnetization as a 
function of the disorder width. 

A difference in the location of the critical point for 
equilibrium and non-equilibrium behavior of the same 
model may appear rather peculiar and one could be 
tempted to ascribe it to finite size corrections. In order 
to clarify this issue, we have solved exactly the model on 
the Bcthe lattice and compared the results for GS and 
DS. While the exponents, as expected, are the same, co- 
inciding with the results of mean-field theory, the critical 
disorder differs in the two cases. Namely the transition 
in the DS occurs at a lower disorder value. Thus there is 
an intermediate parameter region where the GS is ferro- 
magnetic but the DS is paramagnetic. In conclusion, the 
solution on the Bcthe lattice corroborates the picture ob- 
tained from simulations in d — 3. From the optimization 
viewpoint, the d = 3 case shows an intermediate phase of 
"bad" correspondence between the GS and DS, exactly 
as in d = 1. This however stops as the Ri DS ^ is ap- 
proached: naturally if both the states are ferromagnetic 
the optimization of the DS is much easier. To further ex- 
plore the question of universality of the two transitions 
in the GS and in the DS, we have computed the distribu- 
tion of the magnetization at the respective critical point, 
R,i DS ^ and Ri GS ^ for different lattice sizes. The distribu- 



tions can again all be collapsed into the same curve. 

Finally, we consider the question of when is it actually 
possible to reach the exact GS via demagnetization. To 
this end, we consider a reverse field history (RFH) algo- 
rithm that allows in principle to construct a field history 
to get to the GS, if possible at all. Studies of the d = 1 
case illuminate the difficulty of optimizing since it turns 
out that for anything but very strong disorders R the 
probability to reach the GS rapidly decays to zero. 

Our main conclusion is that, in general, demagnetiza- 
tion is not a good technique for reaching states that are 
truly close to the equilibrium, except in cases where the 
outcome is clearly similar from the very beginning (FM 
states and PM states where the disorder is strong) . This 
holds for both the energy of the states and also for the 
spin configurations. A simple formulation is that, since 
the DS is not optimized well in terms of the locations 
of domain walls, it has an excess random field (Zeeman) 
energy. 

The paper is organized as follows: in section II we de- 
fine the model and discuss its numerical treatment. In 
sec. Ill we analyze the one-dimensional case, analytically 
and numerically. Section IV is devoted to the behavior 
around the disorder induced transition in d = 3 and on 
the Bethe lattice. Section V demonstrates the RFH al- 
gorithm, together with numerical studies. Conclusions 
are reported in section VI. An account of some of these 
results was briefly reported in Ref. [34]. 

II. THE RANDOM-FIELD ISING MODEL 

In the RFIM, a spin Sj = ±1 is assigned to each 
site i of a d— dimensional lattice. The spins are cou- 
pled to their nearest-neighbors spins by a ferromagnetic 
interaction of strength J and to the external field H. 
In addition, to each site of the lattice it is associated 
a random field hi taken from a Gaussian probability 
p(h) = cxp(-h 2 /2R 2 )/V^R, with variance R. The 
Hamiltonian thus reads 

H = -Y,J*i*3-Y,{H + hi)8i, (1) 

where the first sum is restricted to nearest-neighbors 
pairs. 

In this paper we will consider only the case of zero 
temperature, both in equilibrium and out of equilibrium. 
The T = equilibrium problem amounts to find the mini- 
mum of 7i for a given realization of the random- fields (i.e. 
the GS) and then eventually perform the thermodynamic 
limit. This problem has been solved exactly in a num- 
ber of simple cases, namely in d — 1 and on the Bethe 
lattice, for particular disorder distributions and studied 
numerically in generic dimensions. 

The RFIM GS is solvable in a polynomial CPU- 
time, with exact combinatorial algorithms. For the one- 
dimensional case, the solution can be found via a map- 
ping to a "shortest path problem" [35] which effectively 



places the domain walls in optimal positions, correspond- 
ing to the global minimum of H. For higher dimensions, 
one starts by noticing that finding the RFIM GS is equiv- 
alent to the min-cut/max-flow problem of combinatorial 
optimization. This can be solved in a variety of ways. 
We use a so-called push-relabel variant of the preflow al- 
gorithm [36]. Such methods, properly implemented, are 
in general slightly sub-linear in their performance as a 
function of the number of spins in the problem. 

For the out of equilibrium case, we need to specify an 
appropriate dynamics, ruling the evolution of the spins. 
We will consider the dynamics proposed in Ref. [37] and 
used in Refs. [13-15] to study the hysteresis loop. At 
each time step the spins align with the local field 

a i =mga(j'^2a j + hi + H), (2) 

3 

until a metastable state is reached. This dynamics can be 
used to obtain the hysteresis loop. The system is started 
from a state with all the spin down Sj = — 1 and then H 
is ramped slowly from H — > — oo to H — ► oo. The limit 
of dHj dt — > can be conveniently obtained by increasing 
the field precisely of the amount necessary to flip the first 
unstable spin. A single spin flip increases the local field 
of the nearest neighboring spins, generating an avalanche 
of flippings. When the systems finds another metastable 
state, the field is increased again. This dynamics obeys 
return-point memory [13]: if the field is increased adi- 
abatically the magnetization only depends on the state 
in which the field was last reversed. This property has 
been exploited in d = 1 [21, 24] and in the Bethe lattice 
[22, 27] to obtain exactly the saturation cycle and the 
minor loops. 

The main hysteresis loop selects a series of metastable 
states, which in principle arc not particularly close to 
the ground state. To obtain low energy states, we per- 
form a demagnetization procedure: the external field 
is changed through a nested succession H = H —* 

Hi -» H 2 -» H n ... -» 0, with H 2n > H 2n+2 > 0, 

H 2n -i < H 2n +\ < and dH = H 2n -H 2n+2 — > 0. A per- 
fect demagnetization can be performed numerically using 
the prescription discussed above to obtain dH/dt — > 0. 
Such a perfect demagnetization is quite expensive compu- 
tationally and it is convenient to perform an approximate 
demagnetization using dH = 1CV 3 . A comparison of the 
states obtained under approximate and perfect demag- 
netization shows negligible differences. 



III. GROUND STATE AND DEMAGNETIZED 
STATE IN ONE DIMENSION 

A. Exact results: ground state 

The GS energy can be computed exactly in d = 1 using 
transfer matrix methods [9] The free energy of a chain of 



length N is given by 



F N = ~ ]n{Z N ) = ~\n(Z+ - Z N ) 



-ln(Z+Z N ) 
(3) 

where Zn is the partition function with free boundary 
conditions and Z^ are the partition functions with the 
spin at site 1 fixed up (down). These functions satisfy the 
following recursive relation: 

Z±=e±^(Z+_ 1 e^ + Z^_ 1 e^) (4) 
The last step in eq.(3) uses the approximation Z^+Z^ ~ 



Z^Z N which holds in the large N limit since Z^ both 

diverge with the ratio Z^/Z^ being finite. From Eq. (4) 
it follows 

Z+Z N = Z+_ 1 Z N _ 1 (2cos\v{(3J) + 2coMWxn)) (5) 



where xn = \n(Z~^/Z N ), which gives for the total free 
energy 



2(5 



ln(2 cosh(/3 J) + 2 cosh(2/3xjv)) (6) 



where xn — \n(Z^/Z N ), so that one can define a free 
energy per site 



/ 



1 



ln(2cosh/3J + 2cosh(2/3a;jv)). 



(7) 



xjv is a stochastic quantity satisfying the equation 

xn = hi +g(x N -i) (8) 

where g(x) = lln((e*+ J ) + l)/(e*)+eW>)) 
When R — > Eq. (8) has a fixed point solution of 
x oo — g{xoo)- It is easy to check that x^ = is the 
only solution for any J and (3 finite, corresponding to 
the absence of a phase transition. 

When R is non-zero xn is a random variable with an 
associated distribution Wm{x), where 

Wn{x)<Ix = Prob(x < xn < x + dx). (9) 

Wn{x) satisfies the recursive functional equation 

W N+1 (x) = f^dhP^x 

f^cLaWNfaWx-h-H-gixi)) (10) 

so that in the thermodynamic limit W x is given by the 
fixed point equation 



W c 



/OG 
cbuWooixJPix-h-H-gix!)). (11) 
-oo 



Once Woo is known, any thermodynamic quantity can be 
computed. In particular, the free energy per spin, which 
is given by 

1 f°° 

(/) =--3 dxW x (x) (cosh(2/3) + cosh 2/3x) . (12) 

P J — oo 



The magnetization at a site of an infinite lattice, is 
given by 

/c \ — zT -z l — 
\ s o; - zi+zi - 



where Z^* are respectively the partition functions with 
the spin at fixed up (down). These arc given by 

z n = e ±0h o{e ± j z + + eT ej z - ){e ±PJ z + + e ^ J zn 

(14) 

where Z^ are the partition functions for the semi-infinite 
right(left) lattice, with the spin at site 1 (—1) fixed 
up(down). This gives 

(s ) = tanh(/?(/i + g(x r ) + g(x t ))). (15) 

Finally, The magnetization for the infinite lattice is ob- 
tained averaging over the quenched variables x r y. 

m = I-oc dhP{h) dx r W N {x r ) 
dxiW N {xi) tanh (J3(h + g(x r ) + g(xi))) . (16) 



B. Exact results: Demagnetized state 

In d = 1 the magnetization and the energy per spin 
as a function of the external field can be derived explic- 
itly through a probabilistic reasoning. We show how to 
get these results on the saturation loop, focusing on the 
lower branch. (The results on the upper branch can be 
obtained by symmetry considerations.) Similar but much 
more involved reasoning can be repeated for any minor 
loop. 

The central quantity to consider, in order to solve for 
the magnetization as a function of the external field H 
on the hysteresis loop, is the conditional probability for a 
spin to be up, conditioned to one of its nearest neighbors 
being down. To calculate this quantity, one can reason 
as follows: fix the spin at site i — 1 down. Define p m (H) 
as the probability for a spin to be up, given that exactly 
m (to = 0, 1, 2) of its neighbors arc up: 



Pm (H) = P(hf > 0) 



dhp(h) , (17) 



where z is the coordination number (z = 2 in d = 1). Fix 
for a moment the spin at site i down as well and look at 
the spin at site i + It will be up with probability Uq 
and down with probability 1 — Uq. The spin at site i will 
flip up with probability p\ when the spin at i + 1 is up 
and pq when it is down. Ultimately, the spin at i will 
be up (conditioned to the spin at i— 1 being down) with 
probability J7o = UqPi + (1 — E/o)Po- It follows 



U 



Po 

1-pi+Po 



(18) 



Once Uo is known, a similar reasoning leads to the (un- 
conditioned) probability p(H) for a spin to be up: Fix 
the spin at site i down. The spin at site i — 1 will be up 
with probability Uo and down with probability I — Uo. 
The same holds for the spin at site i+1. Thus 

p(H) = U 2 p 2 + 2U (1 - U ) Pl + (1 - U ) 2 po, (19) 

from which the magnetization is obtained as m(H) = 
2p(H) - 1. 
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FIG. 1: The energy of the GS is compared with the one of 
DS. The values are computed exactly in d = 1 as a function 
of the disorder width R. 
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FIG. 2: The energy difference between the GS and the DS 
computed exactly in d = 1 as a function of the disorder width 
R. 

The energy per spin on the saturation loop is obtained 
as follows. Due to translational invariance: 

E =^~ = -J{siSi+i) - H(8i) - (h iSi ). (20) 



To calculate the spin-spin correlation (sjSj+i) we intro- 
duce the probabilities <&++, $ + ~, $>-+, $ for adja- 
cent spins to be respectively up-up, up-down, down- 
up, and down-down. These quantities are not indepen- 
dent, since they have to satisfy the obvious identities: 
$+-=$-+, $+++$+- = p(H), and $—+$+- = l-p(H). 
Thus it is sufficient to calculate one of them, for ex- 
ample This is done by separating the four con- 
tributions from the possible boundary conditions deter- 
mined by the values of the spins at sites i — 1 and i + 2: 
When they are both down, the probability for the cou- 
ple of spins at sites i and i + 1 to be both down is 
Ug (1— pi(H)) 2 , when one is up and the other one is down 
it is 2U (l - U )(l - pi (H))(l -p {H)), and when both 
of them are up it is (1 — Uo) 2 (l — po(H)) 2 . Adding up the 
four contributions one gets = (1 — Uo) 2 - This fixes 
the other probabilities to be <f> + ~ = <I>~ + = 2p-l+(l-U ) 2 , 
and <I> ++ = 1—p — (1 — Uo) 2 - Thus, the spin-spin corre- 
lation is 

(siSi+i) = $ ++ + - 2*+-= 4 (p - (1 - U ) 2 ) - 3 . 

(21) 

The average value (ft»s») can be obtained by averaging 
over the held h the product of h times the average value 
of the spin Sj over the local fields other then hi, once the 
field at i is fixed at the value h : 



(hiSi) 



-\-oc 



dh p(h )h (si\h ). 



(22) 



The conditional average (si\ti) is given by 2p(H\h') — l 
where p(H\h ) is the conditional probability for a spin 
to be up at an external field H, given that its local ran- 
dom field is fixed at the value h . From Eq. (19) this is 
trivially given by 



p(H\h) = Ug9(h' + H + 2 J) 

+ 2£/ (l - U )6(h' + H) 
+ (1 - U ) 2 9(h + H - 2 J) , 



(23) 



/ + 0O 
dhp{h)h 
■if — 13 



which finally gives 
(h iSi ) = 2U 2 

+ AU (l-U ) f + dh'p(h')h 

J -H 

/+oo 
dh'p(h')h -h'. (24) 
■H + 2.7 

In particular, for a Gaussian distribution with h' = and 
variance R the integrals can be performed analytically 
and the result is 



(hiSi) 



Re 



2U^e^ cosh (2JH/R 2 ) 



+ e 



^>(l-2U 2 ) + 2Uo{l-Uo) .(25) 



The energy per site on the lower branch of the saturation 
loop is in general given by 

E{H) = -4 J (p(H) - (1 - Uaf) + 3 J - H(2p(H) - 1) 

/ + °° / / / 
dh p(h )h 
-H-2.J 

/+oo 
dh p(h)h 
H 

/+oo 
dh P (h )h - K. 
-Jf + 2.7 



(26) 



Similar but much more involved reasonings can be re- 
peated for any minor loop - eventually for a series of 
nested loops leading to the demagnetized state - pro- 
viding a series of recursive equations for the magneti- 
zation, the spin-spin, and the spin-field correlations, 
which are the quantities needed to compute the en- 
ergy. If the external field is changed through a nested 

succession H = H —* Hi —* H 2 — > H n ... — ► 0, 

with H 2n > H 2n +2 > 0, H2n-i < H 2n +i < and 
dH = H 2n — H 2n+2 — » 0, the spin-spin correlations are 
given recursively by 



(siS i+ i)H 2n - (SiSj+i) iJ 2n _! = 

W* n {p 2 {H 2n )-p 2 {H 2n _i)) 
-4D% n _i(po(H 2n )-p (H 2n _i)) 



(27) 



where C4 and Z?fc are respectively the probabilities for 
a spin to be up(down) conditioned to one of its neigh- 
bors being down, and satisfy in turn a set of recursive 
equations. Similar equations hold for magnetization and 
spin-field correlation, leading to a complicated recursive 
formula for the energy The results of such calculations 
are shown in Figs. (1, 2), where the energy of the demag- 
netized state is compared with the energy of the ground 
state evaluated in the previous section. 



C. Simulations: how optimized is the 
demagnetized state? 

In one dimension the comparison of the DS and the 
GS is the easiest since the domain walls are just point- 
like. For the GS we know that it is optimized such that 
all the large enough local random field fluctuations nu- 
cleate domains of the same sign. The rest of the random 
landscape is split up into regions that align themselves 
with such fluctuations depending on the sign of the ran- 



dom field excess, J2 i 



zGrcgion 



hi. As a result the Zeeman 



energy of domains is linear in domain size, Ez ~ Id, and 
the asymptotic mean domain length follows the Imry-Ma 
prediction (Igs) ~ 1/-R 2 - Moreover since the random 
landscape has a finite correlation length the domain size 
distribution is exponential [35]. 

Any qualitative differences in the DS will follow from 
three separate mechanisms: 1) shifts of domain walls, 
2) creation of domains inside intact GS domains and 3) 



destruction of GS domains (Fig. 3). From the point of 
view of "optimization" the first one is of trivial concern, 
since it would have little effect e.g. on the scaling of 
Hz,ds- The second one is more detrimental if the energy 
difference to the GS is considered. In addition to the 
cost of the two domain walls it subtracts a contribution 
from the Zeeman energy of the domain that persists and 
surrounds (in the GS) the one that is not created in the 
DS. The third one would make the largest change to the 
total energy, since for Id S> 1 the energy of a domain 
consists mostly of its Zeeman energy. 

+ + + - - -- -- + + + groundstate 

( 1 

+ ++;++ j- + + + DW shift 

j I 

+ + + - -;++!- - + + + nucleated 
.ZZZZZ.. droplet 

+ + + + + + + ;+ + + destroyed 

■ i domain 

FIG. 3: An illustration of the possible mechanisms for the 
deviations between GS and DS. 
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FIG. 4: The average change in the spin-spin overlap between 
the GS and the DS (Aq) and the contribution to that from 
completely "destroyed" GS domains (Aq destr ), as a function 
of R. 

Numerical studies of the DS domain structure indicate 
that with decreasing R the average domain size increases 
faster than in the GS, while the size distribution P(ld) 
remains exponential. This is accompanied by a reduction 
in the overlap q = ({o'gsO'ds} + l)/2 between these two 
states. For R large the overlap is close to unity; strong 
local fields hi align the spins in the same way regardless 
of the mechanism by which the spin state is created. For 
R small the local field is no longer strongly correlated 
with the orientation of the spin, and thus whether the 
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FIG. 5: The average overlap of a DS domain of size I with 
the GS domain spin state at the same locations for R = 0.5, 
0.6, 0.8, 0.9, 1.0. The overlap decreases with ldomain,DS- 
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FIG. 6: The Zeeman energy of DS domains as a function of 
the size, ldomain,DS- The black circles mark the average DS 
domain size for a given R. The two lines above and below 
the data indicate optimal, linear (GS-like) scaling and the 
Imry-Ma-like /^-scaling, respectively. 



GS and DS are locally aligned depends on how optimized 
the latter is. 

The fundamental mechanism for the deviations be- 
tween the states seems for R small to be the "destruction" 
of GS domains (see Fig. 3 again). This is demonstrated 
in Fig. 4 by depicting the change Aq in the overlap that 
comes solely from missing GS domains. The conclusion 
from this dominance is that the demagnetized states typi- 
cally miss regions in which the integrated field fluctuation 
is large which as such leads in the GS to the formation of 
GS domain. Therefore the overlap should get smaller the 
larger the scale-length on which one compares the DS and 
GS is, and this is confirmed by Fig. 5 which shows the 



overlap between a DS domain and the GS as a function 
of the length of the DS domain. 

The importance of such destroyed domains can also 
be seen in the total contribution to the energy difference 
between the DS and GS. For R small this is again domi- 
nated by the missing GS domains. In general the differ- 
ence between the energies of the GS and DS derives from 
the combination of domain walls and Zeeman energy. 
Fig. 6 shows that for Id small the DS domains do not 
have much Zeeman energy. This changes if Id is larger, 
in which regime the scaling approaches the Imry-Ma -like 
scaling (/"' 5 ). The implication is that the field energy of 
large domains in the DS self-averages, and comes from a 
sum of random contributions (ie. the domains contain re- 
gions where the actual random field sum is opposite to the 
spin orientation, such as the missing GS domains). The 
cross-over between the small /^-behavior and the asymp- 
totic scaling is located close to (ld)DS- 

IV. AROUND THE DISORDER INDUCED 
TRANSITION 

A. Simulations in d — 3 

The RFIM displays a disorder induced phase transition 
both in the GS and in the hysteresis loop, which can also 
be observed analyzing the DS [21, 22, 39]. If the GS and 
the DS are always paramagnetic, the transition is absent 
and thus we perform numerical simulations in d = 3. Our 
aim is to characterize the difference between DS and GS 
around the disorder induced transition. 

In d = 3 for low disorder, the GS is ferromagnetic, 
while for higher disorder it becomes paramagnetic. For 
Gaussian disorder, the transition point has been located 

(GS) 

numerically at Ri ~ 2.28. It is possible to define the 
usual set of critical exponents characterizing the phase 
transition and compute the values by exact GS calcula- 
tions. For instance, the magnetization M = (\m\), with 
m = J2i S i/N, scales close to the transition point as 

M = Ar p , (28) 

where r = (R — R c )/R c is the reduced order parame- 
ter and A is a non-universal constant. The correlation 
length defines another exponent £ = (Br)~ u -where B 
is another non-universal constant- which rules the finite 
size scaling of the model 

M = AL-P' u f (BL^^R - R c )/R c ) . (29) 

Simulations yield v (GS ^ ~ 1.17 and / 3 (GS) = 0.02. 

A disorder induced transition is also found in the hys- 
teresis loop. At low disorder the loop shows a macro- 
scopic jump, which disappears at a critical value for the 
disorder. This transition reflects itself in the DS, which 
is ferromagnetic when the main loop has a jump and is 
paramagnetic otherwise. The transition point has been 



obtained numerically in d = 3 as Ri° S ^ ~ 2.16 and the 
critical exponents have been measured. In particular, 
Ref. [39] reports data collapses with @(ds) = 0-04 and 
y {DS) = 1-41. While there is strong evidence that the 
exponents measured in the DS should be equal to those 
measured on the main loop, the relation with the equi- 
librium transition is not clear. 

We notice first that numerical simulations reported 
in the literature indicate that the transition appears at 
slightly different locations in the GS and in the DS. Hart- 
mann and Nowak report R GS = 2.29 ± 0.04 for the GS 
with system size up to L = 80, Hartmann and Young 

f GS) 

refine this value to R c = 2.28 ± 0.01 with sizes up to 
L = 96, which is also confirmed by Middleton and Fisher 
which estimate R { C GS) = 2.27 ± 0.04. For the hysteresis 
loop the best estimate is R c = 2.16 ± 0.03, with sys- 
tem sizes up to L — 320 and a similar value for the DS 
[21, 39]. Thus, unless strong finite size effects take place, 
one is tempted to conclude that the two transitions take 
place at two different values of R c . 
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FIG. 7: The magnetization of the GS and the DS in d = 3 for 
different system sizes L and disorder R. 
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FIG. 8: Numerical results in d = 3: The magnetization can be 
collapsed using R c = 2.28 (GS) and R c = 2.16 (DS), l/v = 
0.73 and (3 = 0.03. The scaling curve is the same for DS 
and GS indicating universal behavior. The values for the 
ratios of the non- universal constants are Ads /Acs = 1 and 
B ds /Bgs = 0.68. 
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Here we analyze the problem again by numerical simu- 
lations, computing the GS and the DS numerically, using 
the same disorder realizations for the two cases. Sim- 
ulations are performed for cubic lattices of linear sizes 
L = 10,20,40,60,80 and the results are averaged over 
several realizations of the random fields. The GS is found 
exactly using a min-cut/max-flow algorithm, while de- 
magnetization is performed approximately with the algo- 
rithm discussed in Ref. [21] with dH = 10~ 3 (see section 
II). In both cases, we compute the average magnetiza- 
tion as a function of the disorder width (see Fig. 7). In 
Fig. 8 we collapse the two sets of data into a single curve, 



using two different values for R c (i.e. R c = 2.2<5 
r( ds ^ = 2.16) but the same values for the exponents (i.e. 
1/v = 0.73 and (5 = 0.03). The best value for the ratio of 
the non- universal constant is found to be Ads /Acs — 1 



and 



FIG. 9: The overlap between the GS and the DS in d = 3 for 
different system sizes. 



and B D s/B G s = 0.68 ± 0.02. The fact that the scal- 
ing function is the same for the two cases is a strong 
indication for universality, going beyond the simple nu- 
merical similarity of the exponents. There is always the 
possibility that in the limit L — > oo Ri GS ^ = Ri DS \ At 
the present stage this hypothesis is not supported by the 
data, since we were not able to collapse all the data into 
a single curve using the same R c . 

Next, we compare the statistical properties of the GS 
and the DS around the transitions. In Fig. 9 we report 
the value of the overlap as a function of R for different 
system sizes. When the disorder is decreased from the 
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FIG. 10: The difference in magnetization between the GS and 
the DS in d = 3 for different system sizes. 



paramagnetic region, the overlap decreases as for d = 1. 
However for low disorder the overlap rapidly increases 
and reaches 1 in the ferromagnetic state. The minimum 
of the overlap is located in the parameter region cor- 
responding to the transitions (i.e. R <~ 2.2 — 2.3). A 
decrease in the overlap around the transition can be ex- 
pected, since for r{ DS ^ < R < r{ GS ^ the GS is ferro- 
magnetic (M > 0) and the DS is paramagnetic (M = 0) 
as it is also apparent plotting the difference in the mag- 
netization (see Fig. 10). 

In summary, three dimensional simulations indicate 
that the transitions in the GS and DS are universal, but 
the critical parameter seems to differ. Consequently the 
GS and DS differ mostly around the transition, while the 
difference is smaller in the paramagnetic and ferromag- 
netic phases. 



The Bethe lattice 



lattice. Then, following the d = 1 case, one can write 



where 



F n _ 1 (j)- — ln2(cosh((3J) + cosh(2#r„(j))) 

(31) 



x n {i) = ^Hz+{i)/z-{i)), 



(32) 



so that the contribution at the free energy from site i is 

f(i) = -— ln(2cosh/3J + 2cosh(2/3a;„(i))). (33) 
2p 



x n (i) is a stochastic quantity satisfying the equation 



x n {i) = hi + ^2 9{x n -i{j)) 
jei(i) 



(34) 



When R — > Eq. (34) has a fixed point solution of 
= (z — l)g(x oa ). = is a solution for any J 

and p. For (3 < (i c = \ In there are also two stable 
solutions ±.Too 7^ corresponding to the appearance of a 
ferromagnetic phase. 
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The RFIM can be solved exactly in the Bethe lattice, 
displaying a disorder induced transition in the GS and in 
the DS [22]. It is thus an interesting case to compare the 
two states around the respective transition directly in the 
thermodynamic limit. We consider here a Bethe lattice 
with coordination z and obtain the GS generalizing the 
d = 1 case as in Ref. [38]. In this case N refers to 
the generation of the lattice, and are the partition 
functions of a branch of generation n with a fixed up 
(down) spin at the central site. The recursion relation 
for the Z± is 

Ztii) = e ±phi II {ZUU)^ J + Z-^U)e^ J ) (30) 

j'€/(») 

where for any given site i the sum over j runs over the 
z — 1 nearest neighbors of i away from the center of the 



FIG. 11: The magnetization of the GS and the DS computed 
exactly on the Bethe lattice with z — 4 in the thermodynamic 
limit, showing the ordering of the critical point (see inset). 
When the data are plotted against the reduced parameter 
(R c — R)/R c the curves superimpose. The result implies that 
for the Bethe lattice Aas = Ads- 

To perform quenched averages one has to solve for the 
probability distribution of W n (x n ), where W n {x)dx = 
Prob(x < x n < x + dx), which satisfies the recursive 
functional equation 

/oo r X 

dhP(h) / dx 1 W n (x 1 ) ■ ■ ■ 
-oo J —oo 

■ ■ ■ / dx z ^W n (x z ^)S(x -h-H -$2g{x k )) ,(35) 



so that in the thermodynamic limit Woo is given by the 
fixed point equation 



point, Ri DS ^ and Ri GS ^ for different lattice sizes N. 



The 



IF 



/OO 
dx{Woo{xx) • • • 
-oo 

/oo f 
-°° fe=i, 

Once W^oo is known, any thermodynamic quantity can 
be computed. In particular, the free energy per spin is 
given again by (12) and the magnetization at the central 
site of an infinite lattice, is given by Eq. (13) where 
are respectively the partition function with the spin at 
fixed up (down). These are given by 

Z U = e ±Ph -Q {e ±^ Z + + e^ J Z^) (37) 

k=l,z 



and Z± 



for k = 1 , • • • z are the partition functions of the 
z branches attached to the central site 0, with the bound- 
ary spin fixed up(down). This gives for the magnetization 
at the central site (s ) 



(s ) = tanh(/3(/i + 9{x k ))) 



(38) 



fc=l,2 



The magnetization for the infinite lattice can then be 
obtained averaging over the quenched variables x r y. 



M 



-I 

/OO 
-oo 



oo J —oo 

dx z W N (x z )t<mh(f3(h + ^ g(x k ))).(39) 

k=l,z 



For a Gaussian random-field distribution the fixed point 
equation can not be solved explicitly and we thus re- 
sort to a numerical integration. We obtain Woo(x) for 
z = 4, and for different values of R, and compute the 
magnetization using Eq. (39). In Fig. 11 we compare 
the magnetization of the GS with the one of the rem- 
nant magnetization in the DS, computed in Ref. [22]. As 
observed in the simulations in d = 3, the transition oc- 
curs at two different locations (see the inset of Fig. 11), 
for z = 4 R ( C DS) = 1.781258... [22] and i?i GS) ~ 1.8375, 
with the mean-field exponent ((3 = 1/2). When plotted 
against (R — R c )/R c the two curves superimpose close 
to the critical point. This indicates that, though not re- 
quired by universality, in the Bethe lattice Acs — ^ds, 
as also found in d = 3. 

To investigate possible finite size scaling we have per- 
formed numerical simulations in the Bethe lattice, fol- 
lowing the method of Ref. [25] . Collapsing the order pa- 
rameter curve as in d = 3, using a scaling form similar 
to Eq. (29), does not appear to be possible in the Bethe 
lattice, because the scaling region is very narrow. Thus 
to test finite size scaling, we have computed the distri- 
bution of the magnetization m at the respective critical 



distributions can all be collapsed into the same curve (see 
Fig. 12), using the form P(\m\) = f(\m\/M)/M. 
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FIG. 12: The distributions of the magnetization in the DS 
and the GS at their respective critical points on the Bethe 
lattice, obtained numerically for different lattice sizes N, can 
be all collapsed together. 



V. REACHING THE GROUND STATE BY 
NON-EQUILIBRIUM DYNAMICS 

After having shown that the GS and the DS corre- 
spond to different microscopic configurations, we investi- 
gate now if the GS spin configuration may be reached by 
a field history other then the ac-demagnetization. The 
answer to this question requires a clarification on the 
relation existing between locally stable states (given by 
the solutions of Eq.(2)), and the spin configurations vis- 
ited along the non-equilibrium dynamics induced by the 
varying field. In fact, not all stable configurations may 
be reached by a field history from saturation. The prob- 
lem has been treated in [40] where it has been shown 
that, given a spin configuration obtained by a field his- 
tory, supposed unknown, the sequence of reversal fields 
that applied to saturation gives back the original state 
can be recovered. For spin systems this inverse function 
is given by an algorithm which is able to construct the 
reverse field history (RFH) [40]. This method is applied 
then to investigate if a given spin configuration may be 
reached by field history from saturation: if a field his- 
tory leading to the state exists the algorithm produce a 
sequence of reversal fields; if no field history exists the 
algorithm enters a recursive loop. The investigation of 
the properties of unreachable states has been recently 
performed and leads to a classification of stable states 
on oriented graphs [41]. The study is performed here for 
the GS spin configuration that, for the RFIM at finite 
size and for a given disorder realization, can be indepen- 



dently derived by exact combinatorial algorithms (as the 
max-flow min-cut). 

A. RFH Algorithm 

Consider the final spin configuration s (the set of N 
Ising spins) resulting after the application of a field his- 
tory ending at H = and consisting in a sequence of 
reversal fields {H} = {Hi, . . . , H n } from the saturated 
state and let us define the function s = f({H}). The set 
of all states obtained this way is defined as the hysteresis 
states (H-states). Due to adiabatic dynamical response 
and return point memory, this state s will contain the 
memory of a subset of the reversal fields. In fact not all 
the reversal fields determine the final state s. For ex- 
ample, in terms of average magnetization, the reversal 
fields which give rise to closed minor loops do not influ- 
ence the final state, i.e. their memory is erased, while 
the memory of the set of reversal fields {Hs} which are 
not erased is contained in the final state. The inverse 
function {Hs} = g(s) allows, starting from a spin config- 
uration s belonging to the H-states ensemble, to obtain 
the set of reversal fields {Hs} which have been actually 
stored in the state and that - if applied as a field history 
- will reproduce the original state, i.e. s = f(g(s)). We 
define this set of reversal fields {Hs} as minimal field 
history. 

The RFH algorithm takes as input a configuration s at 
H = and gives as output - when it exists - the reversal 
field history from saturation to the state s. The formu- 
lation of the algorithm is based on the order-preserving 
character of the dynamics [13], and is therefore, appli- 
cable to a wide range of models beyond the RFIM. An 
interesting result of the RFH algorithm is obtained when 
it is applied to a state s not belonging to the H-states 
(i.e. where no field history exists). The iterated search 
for the reversal field sequence enters an iteration and, in 
this case, it can be shown that no field history leading to 
the state exists. 



B. Simulation results in Id 

The RFH algorithm was applied to explore the possi- 
bility to reach the GS by non equilibrium dynamics by 
the numerical study of the RFIM in one dimension with 
periodic boundary conditions. We performed our investi- 
gations on systems having N = 5000 spins, averaging the 
results for 100 different realizations of the same disorder 
R. The GS was obtained by the max-flow min-cut pro- 
cedure for each realization and the RFH algorithm was 
applied. At each disorder value R, the fraction fas of 
the realizations in which the GS resulted to be reachable 
was computed. For comparison the same procedure was 
applied starting from locally stable states generated by 
random sampling the set of local minima. The results 
are shown in Fig. 13. 



As a first finding the GS does not result to be sys- 
tematically field reachable and the fraction depends on 
the disorder. One may conclude that the fact that the 
GS is sometimes reachable is a pure effect of the finite 
system size. However, also for the random states the 
fraction of found states ^rnd sensibly changes with R, 
but following a different curve. If there was no correla- 
tion between GS and the H-states the two curves would 
be coincident. The dependence of iRND on R reflects the 
fact that the number of H-states depends on the disorder 
value and the system size [42], and only at large disorder, 
where the number of locally stable states decreases, the 
ratio between H-states and stable states is significantly 
greater then zero. 




FIG. 13: Fraction of reachable states (averages over 100 real- 
izations of disorder) diamonds: fraction Jrnd ; circles: frac- 
tion fas 



VI. CONCLUSIONS 

For disordered systems like the random field Ising 
model one would be interested in both universality in sta- 
tistical properties and in the question how to "optimize" 
in the case of a sample with a given distribution of the 
impurities. In this paper we have studied this problem 
in detail, by comparing the demagnetized and ground 
states. Our main findings are the following: First, the 
character of the GS is such that it is globally optimized, 
and the demagnetization procedure does not perform well 
unless the optimization problem is rather trivial. This is 
slightly surprising since the conclusion holds in particu- 
lar if the RFIM GS is paramagnetic. Then the DS does 
not manage to find the right spin configuration, so that 
as seen in the d = 1 case many of the domains of the GS 
do not appear in the DS. 

Second, in d = 3 (and with the aid of the Bcthe 



lattice solution), it can be demonstrated that the ex- 
istence of a phase transition for both the DS and GS 
makes the "phase diagram" of optimization to show a 
regime where the outcome is less optimal: in the para- 
magnetic phase of the DS, where the GS is already fer- 
romagnetic since the critical thresholds are ordered such 
that R ( c GS) > R ( c DS) = 1.84. In this regime DS and 
GS are expected to differ strongly in the thermodynamic 
limit. We also provide numerical evidence that the d = 3 
transition appears to have the same critical exponents in 
both the GS and DS [43]. This can be considered both 
surprising - there being no exact field theoretical way 
of treating the d = 3 phase transition - and expected, 
since the functional renormalization calculations in spite 
of their shortcomings indicate that the actions are the 
same [14]. It seems intriguing that such universality is 
met exactly in the limit where the "optimized" character 
of the DS changes. 

The results indicate that for the particular system at 
hand, where the disorder couples directly to the expected 
magnetization, "local" optimization methods have diffi- 
culties. Of course, as in "hysteretic optimization", one 
can perturb or "shake" the state obtained from the DS 
procedure to try to still lower the energy These attempts 
are of course usually just heuristic. In the case of the 



RFIM, the joint approach of optimizing by the DS and 
computing the GS exactly allows to understand better 
similarities and differences between equilibrium and low 
energy non-equilibrium states. 

In addition to the ferromagnetic RFIM model, one can 
consider other systems where two disorder induced phase 
transitions exist. Numerical simulations and analytical 
results have shown that a disorder induced transition in 
the hysteresis loop can be observed in the random bond 
Ising model [44], in the random field O(N) model [45], in 
the random anisotropy model [46, 47] and in the random 
Blumc-Emery-Griffith model [44]. All these systems dis- 
play as well a transition in equilibrium and it would be 
interesting to compare their DS and GS. 

Interfaces in quenched disorder would provide another 
interesting example, since the roughness exponent typi- 
cally differs in and out of equilibrium (i.e. at the depin- 
ning threshold) [4]. It would be interesting to measure 
the roughness of an interface after a demagnetization cy- 
cle (i.e. after the field driving the interface is cycled with 
decreasing amplitude), and compare its properties with 
those of the ground state interface. Finally, there is the 
issue of energetics of excitations in the respective ensem- 
bles: the universality of exponents and scaling functions 
would seem to imply that these also scale similarly. 
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